!>\brief
!! Simple sinusoidal flow in a box, useful to test the Stokes code
!! with a nonpolinomial solution.
!!
!! \n
!!
!! The computational domain is the box \f$\prod_{j=1}^d [0,L_j]\f$,
!! with homogeneous Dirichlet boundary conditions. The viscosity is
!! constant: \f$\nu=1\f$. We then define
!! \f{displaymath}{
!!  \Phi(t;a) = \left(\sin\frac{\pi}{a}t\right)^2, \qquad
!!  \Psi(t;a)=\sin\frac{\pi}{a}t\cos\frac{\pi}{a}t,
!! \f}
!! so that
!! \f{displaymath}{
!!  \partial_{t}\Phi(t;a)=\frac{2\pi}{a}\Psi(t;a), \qquad
!!  \partial_{t}\Psi(t;a)=\frac{\pi}{a}\left(1-2\Phi(t;a)\right)
!! \f}
!! and
!! \f{displaymath}{
!!  \partial_{tt}\Phi(t;a)=-\frac{4\pi^2}{a^2}\Phi(t;a)+\frac{2\pi^2}{a^2}, \qquad
!!  \partial_{tt}\Psi(t;a)=-\frac{4\pi^2}{a^2}\Psi(t;a).
!! \f}
!! The right hand side \f$\underline{f}\f$ is chosen in order to have
!! the analytic solution
!! \f{displaymath}{
!!   u_i = U_0\,\alpha_i\Phi(x_i;L_i)\prod_{j\neq i}\Psi(x_j;L_j), \qquad
!!   p = P_0\, \prod_{j=1}^d\Psi(x_j,L_j),
!! \f}
!! with \f$U_0\f$ being a positive constant, \f$P_0=2\pi U_0d\f$ and
!! \f{displaymath}{
!!  \alpha_i = \left\{\begin{array}{cl}
!!   L_i, & i<d \\\
!!   -(d-1)L_d, & i=d.
!!  \end{array}\right.
!! \f}
!! This yields
!! \f{displaymath}{
!!  f_i = 2\pi^2 U_0\left\{2\left[
!!  \left(\sum_{j=1}^d\frac{1}{L_j^2}\right)\alpha_i
!!  -\frac{d}{L_i}\right]\Phi(x_i;L_i)
!!  +\frac{1}{L_i}\left(d-\frac{\alpha_i}{L_i}\right)\right\}
!!  \prod_{j\neq i}\Psi(x_j,L_j).
!! \f}
!! Notice that for \f$L_i=const\f$ the forcing term becomes
!! \f{displaymath}{
!!  f_i = \left\{\begin{array}{cl}
!!  \displaystyle
!!  \frac{2\pi^2 U_0}{L} \left(d-1\right) \prod_{j\neq i}\Psi(x_j,L_j), & i<d
!!  \\\
!!  \displaystyle
!!  \frac{2\pi^2 U_0}{L}\left\{-2d^2 \Phi(x_i;L_i)
!!  +\left(2d-1\right)\right\} \prod_{j\neq i}\Psi(x_j,L_j), & i=d.
!!  \end{array}\right.
!! \f}
!<----------------------------------------------------------------------
module mod_boxflow_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_boxflow_test_constructor, &
   mod_boxflow_test_destructor,  &
   mod_boxflow_test_initialized, &
   test_name, &
   test_description,&
   coeff_visc,&
   coeff_f,   &
   coeff_dir, &
   coeff_u,   &
   coeff_gradu,&
   coeff_p,   &
   coeff_gradp

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter ::    &
   test_name = "boxflow"
 ! private members
 real(wp), parameter :: &
   ll = 1.0_wp, & ! size of the domain: L_j = L = const
   u0 = 1.0_wp, &
   pi = 3.1415926535897932384626433832795028_wp

! Module variables

 ! public members
 character(len=100), protected ::    &
   test_description(1)
 logical, protected ::               &
   mod_boxflow_test_initialized = .false.
 ! private members
 integer :: d ! dimension
 real(wp) :: p0
 real(wp), allocatable :: l(:), alpha(:), c1(:), c2(:)
 character(len=*), parameter :: &
   this_mod_name = 'mod_boxflow_test'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_boxflow_test_constructor(id)
  integer, intent(in) :: id

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_boxflow_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   d = id

   allocate(l(d),alpha(d),c1(d),c2(d))
   l = ll
   alpha(1:d-1) = l(1:d-1)
   alpha(  d  ) = -real(d-1,wp)*l(d)

   c1 = 2.0_wp*( sum(1.0_wp/l**2)*alpha - real(d,wp)/l )
   c2 = 1.0_wp/l*( real(d,wp) - alpha/l )
   p0 = 2.0_wp*pi*u0*real(d,wp)

   write(test_description(1),'(a,i1)') &
     'Boxflow test case, dimension d = ',d

   mod_boxflow_test_initialized = .true.
 end subroutine mod_boxflow_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_boxflow_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_boxflow_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(l,alpha,c1,c2)

   mod_boxflow_test_initialized = .false.
 end subroutine mod_boxflow_test_destructor

!-----------------------------------------------------------------------

 pure function coeff_visc(x) result(nu_x)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: nu_x(size(x,2))

   nu_x = 1.0_wp
 end function coeff_visc
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,1),size(x,2))

  integer :: i, j
  real(wp) :: phi(size(x,1),size(x,2)), psi(size(x,1),size(x,2)), &
              p_psi_skip(size(x,1),size(x,2))

   call compute_phi_psi(x,phi,psi)

   p_psi_skip = 1.0_wp
   do i=1,d
     do j=1,d
       if(j.ne.i) p_psi_skip(i,:) = p_psi_skip(i,:) * psi(j,:)
     enddo
   enddo

   do i=1,d
     f(i,:) = 2.0_wp*(pi**2)*u0                             &
               * ( c1(i)*phi(i,:) + c2(i) ) * p_psi_skip(i,:)
   enddo

 end function coeff_f
 
!-----------------------------------------------------------------------

 pure function coeff_dir(x,breg) result(d)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: d(size(x,1),size(x,2))

   d = coeff_u(x)
 end function coeff_dir
 
!-----------------------------------------------------------------------

 pure function coeff_u(x) result(u)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: u(size(x,1),size(x,2))

  integer :: i, j
  real(wp) :: phi(size(x,1),size(x,2)), psi(size(x,1),size(x,2)), &
              p_psi_skip(size(x,1),size(x,2))

   call compute_phi_psi(x,phi,psi)

   p_psi_skip = 1.0_wp
   do i=1,d
     do j=1,d
       if(j.ne.i) p_psi_skip(i,:) = p_psi_skip(i,:) * psi(j,:)
     enddo
   enddo

   do i=1,d
     u(i,:) = u0*alpha(i) * phi(i,:) * p_psi_skip(i,:)
   enddo

 end function coeff_u
 
!-----------------------------------------------------------------------

 pure function coeff_gradu(x) result(gu)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: gu(size(x,1),size(x,1),size(x,2))

  integer :: i, j, k
  real(wp) :: phi(size(x,1),size(x,2)), psi(size(x,1),size(x,2)), &
              p_psi_skip(size(x,1),size(x,1),size(x,2))

   call compute_phi_psi(x,phi,psi)

   p_psi_skip = 1.0_wp
   do i=1,d
     do j=1,i-1
       do k=1,d
         if((k.ne.i).and.(k.ne.j)) &
           p_psi_skip(i,j,:) = p_psi_skip(i,j,:) * psi(k,:)
       enddo
       p_psi_skip(j,i,:) = p_psi_skip(i,j,:) ! symmetric
     enddo
   enddo

   do i=1,d
     do j=1,d
       if(j.eq.i) then
         gu(i,j,:) = 2.0_wp*pi*u0*alpha(i)/l(i) * product(psi,1)
       else
         gu(i,j,:) = pi*u0*alpha(i)/l(j)                         &
                    *phi(i,:) * ( 1.0_wp - 2.0_wp*phi(j,:) )     &
                    *p_psi_skip(i,j,:)
       endif
     enddo
   enddo

 end function coeff_gradu
 
!-----------------------------------------------------------------------

 pure function coeff_p(x) result(p)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: p(size(x,2))

  real(wp) :: phi(size(x,1),size(x,2)), psi(size(x,1),size(x,2))

   call compute_phi_psi(x,phi,psi)

   p = p0*product(psi,1)

 end function coeff_p
 
!-----------------------------------------------------------------------

 pure function coeff_gradp(x) result(gp)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: gp(size(x,1),size(x,2))

  integer :: i, j
  real(wp) :: phi(size(x,1),size(x,2)), psi(size(x,1),size(x,2)), &
              p_psi_skip(size(x,1),size(x,2))

   call compute_phi_psi(x,phi,psi)

   p_psi_skip = 1.0_wp
   do i=1,d
     do j=1,d
       if(j.ne.i) p_psi_skip(i,:) = p_psi_skip(i,:) * psi(j,:)
     enddo
   enddo

   do i=1,d
     gp(i,:) = pi*p0/l(i) * (1.0_wp-2.0_wp*phi(i,:)) * p_psi_skip(i,:)
   enddo

 end function coeff_gradp
 
!-----------------------------------------------------------------------
 
 pure subroutine compute_phi_psi(x,phi,psi)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(out) :: phi(:,:), psi(:,:)

   integer :: i
   real(wp) :: xl(size(x,2)), st(size(x,2))
 
   do i=1,d
     xl = pi*x(i,:)/l(i)
     st = sin(xl)

     phi(i,:) = st**2
     psi(i,:) = st*cos(xl)
   enddo
 
 end subroutine compute_phi_psi
 
!-----------------------------------------------------------------------

end module mod_boxflow_test

